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Thermal production of J/tp within quark gluon plasma is reconsidered. We show that if screening 
effects are not strong enough, the "in-plasma born" J/tp's would show up as a peak in the Feynman 
momentum distribution at xf = 0. 
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Over the past fifteen years charmonium suppression has been considered as one of the best signatures of quark 
gluon plasma (QGP) formation. Recently this belief was questioned by some works. J/tp suppression in relativistic 
heavy ion collisions is based on the following simple argument, presented in the original work by Matsui and Satz [1]. 
Lattice simulations show that the heavy quark-antiquark potential becomes screened at high enough temperatures. 
As a consequence, when the range of the potential becomes smaller than the J/tp Bohr radius the bound state can 
no longer be formed. On the other hand, detailed simulations [2] of a population of cc pairs traversing the plasma 
' suggested that, given the large number of such pairs, the recombination effect of the pairs into charmonium Coulomb 
bound states is non-negligible and might even lead to an enhancement of J/tp production. This conclusion received 
■ support from the calculations of [3] , where a two component model for J/tp production was proposed. 

Taking the existing calculations seriously, it is no longer clear that an overall suppression of the number of J/tp's will 
be a signature of QGP. A more detailed analysis is required and a more complex pattern can emerge. In particular, 
we might have suppression in some regions of the phase space and enhancement in others. Indeed, this is the result 
of the analysis presented in [4]. In that work the authors study the fate of cc pairs produced in the early stage of 
heavy ion collisions, comparing the case where they have to traverse QGP with the case where they have to traverse 
ordinary nuclear matter. One of the merits of that paper is to emphasize that very interesting information can be 
q ■ extracted from the scaled momentum (xp) distribution of the produced J/tp's. 

Motivated by these observations, in this work we address the J/tp Feynman momentum distributions. In some 
aspects we follow refs. [3] and [4] with one important difference: we include J/tp production within the plasma. This 
is usually neglected because the number of cc pairs produced in the plasma is believed to be small. However, as it will 
• i-h , be discussed in section III, a closer look into the existing estimates shows discrepancies of two orders of magnitude. 
1 In [5], for example, it was estimated to be of the order of 1% of the total number of charm quark pairs. In [6], with 
the inclusion of thermal parton masses this fraction was estimated to reach 20 — 30%. Finally, in [7] this number 
could, in some cases, be equal to the number of directly produced pairs. No systematic effort was made to reconcile 
these different estimates mainly because the screening mechanism was believed to destroy all bound states. Therefore 
in-plasma charm production has been an issue of open charm physics whereas the main focus has always been on 
hidden charm production and suppression. 

Although well established in several lattice calculations [8] , the color screening mechanism may be not so effective 
in real collisions. The existing lattice results are valid for an infinite mass, static, charge-anticharge pair and we 
know that this is not very close to the real situation, where charges are not so heavy and are moving. Recently, a 
calculation of the screening effect in moving charges [9] showed that, in this case the screening is not so strong. With 
the improvements of lattice gauge determinations [8] of the temperature-dependent quark-antiquark potential new 
calculations of the quarkonium spectrum in QGP were performed. In [10] it was found that at up to temperatures of 
the order of 1.10 T c the J/tp can survive as a bound state in the plasma. At higher temperatures the potential is too 
shallow to hold a bound charmonium state. 

We shall assume that screening is not so strong and, in a narrow temperature window just above the phase transition, 
will allow for charmonium (Coulomb) bound states which "survive" from QGP. More precisely, we will allow up to 
one in hundred bound states to survive in the region where they should be screened and destroyed. As it will be seen, 
these "survivors" tend to escape with xp — 0, thus giving us a new kind of QGP signal. Moreover, this peak at low 
xp will fill the dip predicted in [4]. 
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We will consider J ftp production in nucleus-nucleus collisions by three different mechanisms: i) direct (primordial 
QCD parton fusion) ii) thermal (statistical coalescence at hadronization) and iii) QGP (in-plasma parton fusion). As 
already discussed in [3] i) and ii) are just two extreme cases. In realistic simulations, for a fixed number of cc pairs 
and of J ftp's, these last are gradually destroyed, giving origin to the "regenerated" component ii). In the particular 
calculation presented in [3] , for central collisions and for high enough energies a near to complete replacement of the 
initial by the final thermal J/ip's was found. Component iii) was considered long time ago [11], before the work of 
Matsui and Satz, and then it was left in oblivion. Once produced the charmonium state will suffer interactions with 
the partons in the plasma, or with hadrons in hot and dense medium and in a later stage with comoving hadrons 
[12,13]. 

Before performing detailed calculations one could try to guess the shape of the outcoming J ftp's xf distribution in 
the three mechanisms mentioned above. QCD production has been studied in the past and the naive extrapolation 
of the spectra measured in p — p and p — A collisions would lead us to conclude that the J ftp's produced in i) should 
follow a distribution of the type [14] : 

~ {i-x F y (i) 



where c ~ 5 — 6. As it will be discussed in the next section, perturbative QCD (pQCD) calculations within the 
framework of the color evaporation model (CEM) confirmed this behavior in leading order and next-to-leading order 
[15]. No explicit calculation of this xf distribution was performed with the color octet model (COM), but a recent 
study of asymmetries in open charm production [16] suggests that a similar smooth behavior at x F — might be found 
for the J ftp differential cross section. Of course, this has to be verified. Thermal production follows a Bose-Einstein 
distribution [17], which integrated in p\ gives: 



exp 



1 



(2) 



where xf and yfs are the Feynman momentum and the c.m.s. nucleon-nucleon invariant energy. Comparing (1) with 
(2), we can see that the latter flattens out at low values of xf, because of the mass mj/^, forming a "plateau". A 
spectrum from mechanism iii) is a little more difficult to predict. It involves the convolution of two distributions of 
the type (2) (for the two colliding massless partons) and therefore we expect it to be more steeply falling with xp. 
As it will be seen, within certain approximations, it is given by: 



This behavior was found with the CEM. It remains to be shown that it is not an artifact of the model and that it 
would also be found in the thermal version of the COM. Comparing (1), (2) and (3) we see that, if QGP is formed 
and if the multiplicity of "in-plasma-formed" J ftp's is large enough, we may expect to see an enhancement (or even 
a peak) at the origin of the xf distribution. This signal is interesting because it was not predicted by any other 
production mechanism. However it relies on a very accurate knowledge of the pQCD, thermal and QGP abundancies, 
which is very difficult, as it will be dicussed in detail later. 

At this point one might argue that after being produced by the usual pQCD mechanism, the cc bound state suffers 
interactions either with the plasma or with comovers in a hadron gas. These interactions (called here Final State 
Interactions) will distort the initial xf distribution of the bound state. In Ref. [4] it was shown that FSI (either with 
a plasma or a hadron gas) will suppress the charmonium yield at xf = 0, giving origin to a dip in the central xp 
region. 

We are suggesting that a peak will arise at xf — as a consequence of in-plasma production. This is a signal which 
will emerge from the background composed by initial pQCD production. If we ignore the final state interactions we 
will calculate a background which already contains an enhancement of J ftp production in the low xp region, making 
it harder to detect the signal. In other words, if we treat correctly the FSI this will lead to a background with a dip 
around xf = instead of an enhancement. Therefore, at least for the sake of our argument, neglecting FSI is being 
conservative. This is the correct procedure in order to isolate the new effect, which we are looking for. Nevertheless, at 
the end of this work we will include FSI quantitatively and investigate their effect on the Jf tp momentum distribution. 

The text is organized in the following way: in section II we describe direct production; in sections III and IV 
we address production in the plasma and final thermal production respectively; in section V we discuss final state 
interactions; in section VI we present numerical results and in the last section we make some concluding remarks. 
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II. INITIAL J/ip PRODUCTION IN NUCLEUS-NUCLEUS COLLISIONS 



Nowadays perturbative QCD calculations of cc production can be found in textbooks. Nevertheless, we shall, 
in what follows, give some formulas to introduce the notation, to stress the role played by shadowing and to obtain 
expressions in terms of Xf, the variable of interest. One popular approach to charm production is the color evaporation 
model [15,18,19]. We are aware of the discrepancies between different implementations of perturbative QCD in 
charm production [18]. In particular, there is another model, namely the color octet model [20], which gives a more 
sophisticated description of the hadronization process in terms of non-perturbative matrix elements. While the COM 
seems to give a very good description of high px data (pr > 2 GeV) , especially those from the Fermilab Tevatron, it is 
not so reliable at smaller pr [21] (pr < 2 GeV), where the bulk of charmonium production takes place and where we 
wish to make a comparison between QCD and plasma production. For our purpose of studying the xf distributions, 
as can be inferred from the calculations presented in [22], the COM prediction is compatible with a smooth behavior 
of the cross section da/dxF at the origin xf — 0, much like the CEM. Therefore, for simplicity, we shall use this latter 
model in what follows. Our results will certainly depend on the choice of approach. Inspite of the uncertainties, we 
will assume that, since we are interested in the low xf region, pQCD is enough and nonperturbative effects are small, 
in sharp contrast to what happens at large xp [23]. 

In the CEM, charmonium is defined kinematically as a cc state with mass below the DD threshold. In leading order 
(LO) the cross section is computed with the use of perturbative QCD for the diagrams of the elementary processes 
qq — ► cc and gg — > cc convoluted with the parton densities in the projectile and in the target. Calling xp the fractional 
momentum of the produced pair (with respect to the momentum of a projectile nucleon in cm frame) and \fs the cm 
energy of a nucleon-nucleon collision, the cross section for production of a cc pair with mass to is just given by: 



do-PP- 



dxFdm 2 



dxi dx2 5{x\X2S — to 2 ) S(xp — x\ + x 2 ) H {x\, x 2 ; to 2 ) 



s y/x F + Am 2 /s 



H(x 01 ,x 02 ;m 2 ) ; x 01fi2 = 7; (±x F + \j x 2 F + 4m 2 /sj 



(4) 



where x\ and x 2 are the nucleon momentum fractions carried respectively by partons in the projetile and target. The 
function H(xi,X2;m 2 ), which represents the convolution of the elementary cross sections and parton densities is given 
by: 



H(x 1 ,x 2 ;m 2 ) = f g (x 1 ,m 2 )f g (x 2 ,m 2 )a g g(m 2 ) 

[f q (xi,m 2 ) fq(x 2 ,m 2 ) + fq(x l7 m 2 ) f q (x 2 , to 2 )] a q g(m 2 ) 



q=u,d,s 



(5) 



with the parton densities fi(x,m 2 ) in the nucleon computed at the scale to 2 = xix 2 s. 

The LO elementary cross sections in terms of the pair invariant mass (to) are given by [24] : 



,(to 2 ) 



na 2 (to 2 ) 
3to 2 



Ami m, 
1 + — ^ + 



TO^ 



TO^ 



In 



1 + A 
1 - A 



7 + 



31m; 



TO" 



(6) 



87ra 2 (TO 2 ) 
27to 2 



2to; 



TO^ 



Am 2 



TO 



1/2 



(7) 



where m c is the mass of the c quark. 

The production cross section of the charmed state i (= J/ip, ip' or Xc.j) is then finally obtained by integrating 
the free pair cross section cc over the invariant mass to starting from the production threshold 2to c up to open charm 
production threshold 2tod: 



d(T p P ^j/i> 



dxi 



= F 



J 4m? 



dm" 



, da pp - 



dxpdm 2 



(8) 



where F j/^ is the fraction of a cc which contains the corresponding cc resonance. 

This model describes well the experimentally measured x f distribution of hidden charm both with LO and NLO 
cross sections, provided that F^° is defined as Ff LO multiplied by a theoretical factor n, which is equal to the ratio 
of the NLO and LO cross sections {F^° w 2, 54%) [25]. 
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In what follows we shall use the CEM to study perturbative J/tp production in nucleus- nucleus collisions at RHIC. 
As it is well known, in nuclear collisions and in processes involving small values of x, shadowing plays an important 
role (see, for example, [26]). In this case expression (5) is rewritten as: 

H AB (xi,x 2 ;m 2 ) = AB < f g (x l7 m 2 ) Rf(xi,m 2 ) f g (x 2 ,m 2 ) Rg (x 2 ,m 2 )a gg (m 2 ) 



E 

q—u,d,s 



/ 9 (a;i,m 2 ) Rf(x 1 , m 2 ) f q (x 2 ,m 2 ) R? {x 2 ,m 2 ) 

a qq (m 2 )\ (9) 



+ fq(xi,m ) Rg (xi,m ) f q (x 2 ,m ) R q (x 2 ,m ) 
where 

R A (x,m 2 ) = f ^ X,m ^ (10) 

with f A (x,m 2 ) being the i parton momentum distribution in a nucleon inside the nucleus A. Replacing H by Hab 
in (4) we obtain the cross section for cc production in A — B collisions: 

da AB ^ cS f°° . , 1 



dx F Jiml s ^Jx 2 F + Am 2 /s 

and the J/ip production cross section as: 



f 1 

= K dm 2 = H AB (x i,x 02 ;m 2 ) (11) 

Jimi s \/xi? + Am z s 



kF^° / dm 2 ; = H AB (x Q i,x 02 ;m 2 ) (12) 

dx F J Ami s V X F + 4m 2 /s 



In a central A + A collision the number of J ftp's produced directly is related to the cross section by [27]: 



dxF 7r R 2 ^ dx} 



(13) 



where R A is the radius of nucleus A. 

The expression above would give the final distribution if there were no absorption due to further interactions with 
partons in a QGP or hadrons in a comoving fireball. The effect of these final state interactions will be discussed later. 

Since we are going to compare the total number of initially produced cc pairs and J ftp's with those produced in 
the plasma it is useful to introduce the compact notation 

1 f 1 , da AA ^ c5 
-Jo 



^co = ^ / dXF (14) 



N J '+= /"dx f % (15) 



for the number of cc pairs and 

j-AA^J/ip 
dire 

<QCD-J o dx 

for the number of J ftp's. 

Expressions (14) and (15) are calculated with the nuclear parton distributions which take shadowing into account. 
While this makes calculations more realistic, it also introduces some model dependence in the results. In order to have 
a baseline for comparisons it is useful to introduce the equivalent definitions of Nq CD and Nqq D without shadowing. 
For central A — A collisions these are: 



A7CC 

1S QCD 



f 1 da pp ^ cc 

= T AA (b = 0) / dxp — : (no shadowing) (16) 

Jo dx F 



and 
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P 1 j-1 'DT) — ► CC 

N QCD = F^°T AA {b = Q) / dx F ^- (no shadowing) (17) 

where T^yi is the usual nuclear overlap function. In Table I we present some quantitative results for (14), (15), (16) 
and (17). As it can be seen, these numbers change with different choices for the charm quark mass and for the parton 
distribution. Inclusion of shadowing reduces the number of charm quark pairs and J/tjj's in about 10 %, choosing 
different parton densities may change these numbers in 40 % but what changes the most our results is the adopted 
value for the charm quark mass. This choice may alter the numbers by a factor 4. Our number of cc pairs is smaller 
but still compatible with other estimates presented, for example, in [18]. On the other hand, the number of J/ip's can 
be estimated taking into account the recent PHENIX data [28] onp-p collisions. These estimates [29] indicate that 
Nq£ D ~ 0.1, which is close to the number of J/ip's in Table I. 



III. J/ip PRODUCTION IN THE PLASMA 
A. The number of cc pairs 

As in the previous section we reproduce below (to a great extent) textbook material [30], with the purpose of 
defining our notation. The derivation presented here is also useful because, relaxing the screening hypothesis, we 
calculate the production of charm bound states in the plasma and we integrate the rate on all variables except xp. 
All this requires some straightforward but not very often shown manipulations. 

The computation of the in-plasma cc pair production rate goes back to the late eighties [31], was discussed in 
short papers (for example [5]), included in comprehensive review articles [32] during the nineties and experienced 
improvements due to advances in thermal field theory [6]. 

Assuming that QGP is formed, we then have a gas of quarks and gluons with momenta obeying respectively Fermi- 
Dirac and Bosc-Einstein distributions that can collide producing cc pairs. This is the same mechanism, which in the 
strange sector causes the strangeness enhancement. In the charm case, there will be always a "charm enhancement" 
but, because of screening, the charmed quarks will mostly form open charm. In the present calculation we will let 
part of the thermally produced charm Coulomb bound states escape. Since the temperature is high, there will be a 
significant number of parton-parton collisions at a cms energy high enough to produce charm quark pairs, which may 
form charmonium states. We will now estimate their production rate using the CEM in a thermal environment. 

The charmonium production rate in the reaction gg — > cc, at temperature T is given by [30] : 



dN^ cS 1 1 



£ I d 3 Pl d 3 p 2 f g (E{) f g (E 2 ) a^ g °(m 2 ) v 12 (18) 



dtd 3 x 2 (2tt) 6 ^ 

where g g is the gluon statistical factor (number of colors x number polarization states) v\ 2 is the relative velocity 
between colliding gluons with energies Ex and E 2 and three momenta pi and p" 2 , <r gg ° is the elementary gluon-gluon 
cross section (6) and f g (Ei) the usual thermal distribution function: 

We now introduce the charm pair four momentum p with help of the delta function S 4 [p — (pi +P2)]: 

J Ajgg^cc 11 [' 

~^p— = 2 J2nf 9 ° J dtd * xd *P^P* /^O ^°(™ 2 ) v ^ 6 *\P (Pi + P*)] ( 20 ) 
In the expression above we have: 

Pi = (Ei,p Tl ,p Zl ) ; p 2 = (E 2 ,p T2 ,p Z2 ) ; p=(E,p T ,Pz) (21) 

m 2 = ( Pl + P2 ) 2 ; v X2 = \-f- (22) 

where m is the invariant mass of the pair. We next decompose the delta function into temporal, transverse and 
longitudinal components: 
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d^d7 z = W 5 * J ^/n^/v,f 9 mf 9{ E 2 )^o { ^W 



x 6[m 2 - (pj + p 2 ) 2 } S 2 [p T - (p Tl +Pt 2 )} S[ Pz - ( Pzi +p Z2 )} 

Making use of the identity 

^d 3 Pl = J d^SipfWEi) = J dE t dp Tl dp z AE 2 -pl-pl]6m 



we arrive at: 



dNsa~ 



2 r poo roc roc 

dmUp T d Pz =Wf£ J dtd3x l ^Jjp^Jjp^El-pl -pl MEl 

/'OC /■ OO /'OC 

x/ d£ 2 / dp T2 / dp Z2 S[E 2 2 ^p 2 t2 -p 2 Z2 }9(E 2 )f g (E 1 )f g (E 2 ) 

JO J -oc J -oo 



x a^(m")m 2 5[m 2 - { Pl + p 2 ) z ]5 z [p T - (p Tl + pt 2 )]S[ Pz - (p Zl +p Z2 )] 
Integrating in p Z2 , pr 1 , Pt 2 and E 2 and denning [5 as the angle between pr and pVi , we finally obtain: 
dN"^ cS - r °° /> °° r2n 



/>OC />CG pZTT P /'OC /-OO 

= 17^6^ / dm l d PT d(3 dtd 3 x dE 1 dp Zl f g (Ei) f g (E — E{) 



x *£V 2 ) ^ <5[™ 2 - 2EE 1 + 2p T (E 2 1 - p 2 zi y/ 2 cos(/3) + 2 PzPzi ] 



99 V- / £ 

x 9^) 0(E - E^ 9(E 2 - p 
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4(2tt) 



|ff(/n)| 



0{E 2 -h\) 



\H{h 2 )\ 



6{E 2 -h 2 ) 



eiE^eiE-Ex) 



where £7, H(hi t2 ) and fti^ are given by: 

£=[™ 2 +Pt+P 2 ] 1/2 ; H(h 1>2 ) = 2 Pz 



2pT cos(/3) 
(^i 2 -^ 2 , 2 ) 1/2 



'1,2 



(23) 



(24) 



(25) 



(26) 



(27) 



hi,2 = 



1 



^(2^ -m 2 ) T {PT cos 2 (/3) 



2 P 2 + PtCos 2 (/3) 
x ^(p^cos 2 ^) +p 2 ) - (m 2 - 2EE 1 ) 2 }} 1 / 2 



(28) 



In the first line of (26) we have S[m 2 — 2EE\ + 2p T (E\ — p 2 ^ 1 ^ 2 cos(/3) + 2p z p Zl ]. It is easy to see that taking either 
Ei ~ p Zl or £?i >> p Zl and using the <5 function to perform the integral in p Zl we obtain a factor l/p z . This factor 
survives and, at the end, gives the 1/xf structure (3) mentioned in the introduction. The actual numerical calculation 
involves no approximation and has a similar behavior. Notice that, because of the kinematical cuts in the invariant 
mass integral, we are now restricting the result to the rate to charmonium states. A similar expression can be derived 
for cc production originating from quark-antiquark annihilation: 



dNW~ 
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\H(h 2 



B{E{)e{E-E 1 ) 



(29) 



where g q ^ is the statistical factor for quarks (antiquarks) (number of colors x number of spin states), a-^°(m 2 ) is 
given by (7) and f q ^(Ei) is the usual quark (antiquark) thermal distribution function: 
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fq,q( E i) ~ eEt / T , 1 



(30) 



In order to account for expansion effects we shall assume that the system cools down following Bjorken hydrody- 
namics, in which temperature and proper time are related by [30]: 



T(r) _ /r„\V3 



(?) 



(31) 



where To is the initial temperature and To is the thermalization time, which marks the beginning of the hydrodynamical 
expansion. With the help of this relation we can perform the following change of variables: 



dtcPx = dx±_ rdydr = -3 dV t T 3 T~ 4 dT 



j dtd 3 x = 3 J 



T T*T- 4 V(T)dT 



Finally we introduce the energy scale y/s (the invariant energy of a single nucleon-nucleon collision): 



Pz = x F - 



(32) 



(33) 



Expressions (26) and (29) can thus be rewritten as: 
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The volume of the system evolves in time according to: 



V(t)=V - 



TO 



V(T) = V 



(35) 



(36) 



where V = tt R\ t , Ra = (1.18 A x l z - 0.45) fm and the k factor was introduced explicitly. 

Summing Eqs. (34) and (35) and integrating over xp we obtain the total number of in-plasma produced cc pairs, 
which we call Nq GP . 



B. Comparison with other works 



We shall now compare the number of cc obtained in this work with the number of charmed pairs obtained in other 
works. There are some known papers on the subject [5-7,32,33]. From the reading of these papers, we conclude that 
there are large discrepancies in the numbers and in the way to obtain them. The sources of these discrepancies are: 

a) initial temperature of the plasma, T (ranging from 300 MeV to 550 MeV) 

b) degree of parton equilibration (described by the fugacity factors) 

c) initial volume and/or thermalization time, Vq and To 

d) total energy contained in the fireball 

e) use or not of a k factor (=2) in computing thermal rates 

f) use or not of temperature dependent a s 

g) mass of the charm quark m c (going from 1.2 to 1.5 GeV) 
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h) use or not of thermal masses for gluons and quarks in the reactions g + g — > c + c and q + q — > c + c. 

Depending on the choices that one has to make in dealing with a) — > h) the final value of in-plasma produced cc 
pairs can change by orders of magnitude, going roughly from Nq GP — 0.02 in [5] to up 15 in [32]. In a comprehensive 
analysis this same variation was found in [7] . In [6] , the authors arrive to the conclusion that in-plasma production 
was only a factor two smaller than the initial production. However in that paper the factor n — 2 was not included 
in calculation. If it were, then Nq GP ~ Nq CD . 

The list of uncertainties given above is not unique and could be enlarged. Each of the items implies taking a decision 
or making an assumption. Since the list is already large it seems that a high precision calculation is hopeless. In fact 
the situation is not so bad since our knowledge will increase both because of new experimental results and because 
of lattice results. RHIC data will impose severe constraints on items from a) to d). For example, taking together 
data from all the four collaborations, we may expect to know with sufficient accuracy the rapidity distribution of 
charged particles, from which we can have a good knowledge of the total energy contained in the fireball (item d)). 
The global analysis of data on rapidity, pr spectra, abundancies, eliptic flow and HBT interferometry, will eventually 
rule out several initial conditions used in hydrodynamical models and we will have much less uncertainty in the 
initial temperature of the plasma. On the theoretical side, the perturbative QCD analysis of other processes like, 
for example, J/tp photoproduction or J/ip production in e~e + collisions may significantly reduce the uncertainty in 
the charm quark mass (item g)). Lattice calculations at finite temperature will hopefully reduce the freedom in the 
choices of items f ) and h) . A review of what we may learn in the field in a near future can be found in the summary 
talk presented by Pisarski in the last Quark Matter meeting [34] . 

We have explicitly investigated the effect of changing the mass of the quark c and the differences which arise when 
we use the coupling constant running with ra 2 : 

aJm 2 ) — 5— (37) 

(33-27V F )l nx |2 2 



QCD 



or with T [6]: 



^= (33-2^^ (38) 



where N F = 3, A QC d = 230 MeV and = 80 MeV. 

In Table II we present our results for the number of "in-plasma born" cc pairs for different values of couplings and 
charm quark masses. All the calculations were done with an initial plasma temperature of To = 550 MeV, a s (M 2 ) 
and a s (T) are given by (37) and (38) respectively. The numbers inside parenthesis correspond to the choice k = 1. 
Otherwise the numbers were obtained with n — 2. 

The fourth line of Table II shows that our results may change by a factor 30 depending on the inputs used. Since 
there is no strongly preferred value for the c quark mass, neither for n or for the functional form of a s , we are not 
able to a priori discard any of these choices and our final results will reflect these uncertainties. In some cases a direct 
comparison with other works is possible. For example, in the second column, comparing the numbers in parenthesis, 
we notice that we obtain nearly five times more pairs than in [6] . A similar excess is observed comparing the numbers 
in parenthesis in the fourth column. Comparing our work with Ref. [6] we can see that the pair production mechanism 
is quite similar but the treatment given to the plasma expansion is different. Whereas we have used the standard 
Bjorkcn hydodynamics, in Ref. [6] a new hydodynamical model was introduced. Comparing the details of both 
approaches we concluded that in Ref. [6] the expansion and cooling of the system is much faster than in Bjorken 
hydrodynamics. Consequently the system stops much earlier to create cc pairs and the final yield will be smaller. We 
think that the expansion of the plasma must be included in any serious calculation, but hydrodynamics is in itself a 
very complex subject. A prudent strategy, which we adopt here, is to study the desired effect (charm production) first 
with the most standard and simple hydrodynamical model and then, as a second step, plugg the charm production 
formalism into a state-of the-art hydrodynamical code. We will leave this last step for the future. The comparison of 
the other results in Table II with other works shows that for similar inputs we obtain numbers which are compatible 
with those presented in Refs. [7] and [32]. 



C. The number of J/^'s 

The number of J/ip's produced in the plasma can be obtained from Eqs. (34) and (35). We must change the upper 
limit of integration in m 2 introducing a kinematical cut-off, i.e., making the replacement 00 — > 4m|,. In doing so, 
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we rename the superscripts in Eqs. (34) and (35) to gg — ► J ftp and qq — > J/V'- Moreover we introduce the CEM 
multiplicative factor Fj/,/, and arrive at: 

= /s ^ [ dx F + dx F ] (39) 

and 

The color evaporation model has an "intrinsic efficiency" given by the fractional factor Fji^ above, which was 
fixed [18,15] in the analysis of p — p reactions to be Fj/^ ~ 0.02 = 2%. It is by no means obvious that the same 
value should hold for A — A collisions. For simplicity, we shall assume that Fj/$ is universal and holds even for 
in-plasma production. The "screening factor" f s is thus the only free number introduced here. It gives the probability 
that a J/ip formed inside the plasma survives the passage through the medium. In other words, f s accounts for 
dynamical screening, in which gluons in the plasma destroy the charmonium bound state. It varies from to 1. A 
literal interpretation of [1] or [8] would imply f s — 0. We will however tolerate a small value of f s and examine the 
consequences. In fact, it is impossible to say what a "realistic" value of f s would be. Our choices, based on some 
numerical estimates, are in the following interval: 

10~ 3 < f s < 1(T 2 (41) 

If f s < 10~ 3 we do not observe any visible effect of the in-plasma production. On the other hand, f s = 10~ 2 , which 
can be interpreted as meaning that 1 % of the in-plasma born J/ip's can survive as bound states, can be taken as an 
upper limit, beyond which, rather than taking small fluctuations into account we would be really challenging the well 
established concept of screening. 



IV. FINAL THERMAL J/ip PRODUCTION 



We will also, for completeness, consider the case where the J/^'s produced in the early stage of the collision 
traverse the plasma being destroyed and then "regenerated" [3]. We call them "thermal". As has been discussed 
in the literature [17,35], regenerated J ftp's follow a thermal distribution with the temperature of the quark-hadron 
phase transition T c ~ 170 MeV: 



AA^J/i> 



thermal 

dXF 



TTyfiVoA 



dp T exp 



— {xls/l + pl + m*^) 1 ' 2 



= const exp 



where the constant in front of the integral will be fixed later through the normalization condition: 

r-l 



A, 



//V _ 



/,/>, 



o 



dN AA^.7/4, 
dx F thermal 

dxp 



(42) 



(43) 



V. FINAL STATE INTERACTIONS 



In sections II, III and IV we discussed J/ip production. After being produced these states will suffer interaction 
with the free partons of the plasma and later with the hadronic comovers. In [4], the final state interactions (FSI) were 
incorporated (to the rate analogous to our (13)) by the introduction of the multiplicative suppression factors Sp®^ 
and Spgf, for interactions with comovers and QGP partons respectively. We shall adopt here the same procedure. In 
fact, in the case of the initial QCD production, since we are using the same formalism, for the sake of comparison we 
shall later borrow the final corrected expression from [4]. In the case of (40), the J/ip dissociation by partons within 
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the plasma is already taken into account by the "screening factor" f s . If the charmonium state survives the plasma 
phase, it will interact with hadronic comovers and (40) will be multiplied by the suppression factor [36,4]: 

(44) 

where n co is the comover density, nj is the freeze-out density (= 1.15/m -2 ) and a co ~ 1 — 4 mb is the charmonium- 
hadron cross section [12,13]. The above expression depends on the impact parameter, on the collision energy and on 
time at which the charmonium interacts with a given comover. Moreover, since n co is a function of the rapidity, it will 
depend also on xf- However, restricting the analysis to central collisions and to a fixed (y/s) energy, the suppression 
factor tends to be constant, as concluded in [3,36,37]. For simplicity we shall take here the average of the values 
quoted in these works: 

S^g M = 0.25 (45) 

Final state interactions with comovers will be included in the thermal production differential rate (42) in the same 
way and with the use of (45). 



oCOM _ 
'-'FSI — ex P 



-cr co n co In 



fo 



VI. FEYNMAN MOMENTUM DISTRIBUTIONS 



In Fig. 1 we show, for Au - Au{y/s = 200 GeV) (la) and Pb - Pb(y/s = 17 GcV) (lb) collisions, the results 
obtained with (13) (dash-dotted line), (39) (solid line) and (42) (dashed line). For the QCD distribution, we have 
used the nuclear parton distribution functions (PDF) parametrized as in [38]. In doing so, the proton PDF's were 
taken from [39] (GRV98 LO). We also used m c = 1.2 GeV, m D = 1.87 GcV and k = 2. With these choices the 
total number of produced J/^'s is determined. In computing (39) we have used r = 0.7 fm and T = 550 MeV, for 
Au + Au and r = 1 fin and T = 300 MeV for Pb + Pb collisions. T c was taken to be 170 MeV. 

In each figure all curves are normalized to the same number of produced J/ip's, which is given by (15). In Fig. 1 
we have used the parton distributions of [38], obtaining Nj/^ ~ 0.09 at RHIC and Nj/^ ~ 0.001 at SPS. 

Although the normalization is still artificial (it overestimates the number of J/ip's produced in QGP) it anticipates 
our main claim: QGP production will create a peak at xp — 0. Even if suppressed, it will leave a signal. Moreover, 
this may happen at RHIC and at SPS as well. 

In order to further investigate the differences between the three J/ip production mechanisms we will study the 
relation between the measured number of J/ip's and the naive expectation based on what we know from pp collisions. 
Experimentally this ratio is easily constructed from the measured Xf spectra in A+A and p+p collsions. Theoretically, 
the numerator can be written in terms of our theoretical prejudices. 

As a starting point, we assume that all J/^'s are produced directly and we show in Fig. 2 the following ratio: 
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RAx F ) = ^ ^ T7T = —, 7TT (46) 

where the denominator is given by the product of (8) with the nuclear overlap function TaaQ> = 0) = 29.9 mb^ 1 at 
RHIC and = 30.4 mb -1 at SPS. The numerators were obtained with (12), (13) and the nuclear parton distribution 
functions (PDF) parametrized as in [38] (solid line) and as in [40] (dashed line). In the first case the proton PDF 
were taken from [39] (GRV98 LO) whereas in the second case they were taken from [41] (MRST LO). With this last 
set of parton distributions N J/4 , ~ 0.06 at RHIC and N J/4 , ~ 0.001 at SPS. 

Ri shows essentially the effect of shadowing in the low x region. We can see that if no effect is present, other than 
shadowing, J/ip spectrum produced in AA collisions is just a constant times the corresponding pp spectrum. As it 
can be seen this constant may change a lot with the PDF. This constancy in the region xf < 0.1 will be important 
for the subsequent discussion. This figure shows also that the value of the constant decreases when we go from SPS 
to RHIC collisions. This same behavior was found in [4]. 

In Figure 3 we present the xf distribution of the in-plasma produced J/ip's. The screening factor will be fixed to 
f s = 10~ 2 . The curves are obtained with expression (39) normalized to the largest (solid line) and to the smallest 
(dashed line) number of J/ip's. These normalization constants correspond to the choices 2.7 and 0.16, in the last line 
of Table II. Fig. 3a and 3b correspond to RHIC and SPS collisions respectively. 

If we believe that all production is mainly direct plus some small plasma component, then (46) aquires the form: 
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AA^J/-4, 



Rn{x F ) = ^ (47) 



T AA (b = oy 



dxT. 



where, for the moment, we neglect the final state interactions discussed in the previous section. If the screening in the 
plasma is so strong that f s — 0, then, Rn — Ri ~ const. If, however, there is a small J/ip survival probability, for 
example f s = 1CP 2 , there will be a noticeable change in Ru. This is illustrated in Fig. 4, where we see a pronounced 
deviation from a flat behavior in the low xf region if there is a contribution from the plasma. This interesting 
feature is due to nature of plasma production, which is peaked at xp ~ 0. In Fig. 4 the two upper lines correspond 
to the largest plasma J/ip production {Nqq P = 2.7) and the lower two lines to the smallest plasma contribution 

(Nqqp — 0.16). The scrrcning factor is kept fixed to f s = 10~ 2 . In all cases solid and dashed lines mean the same 
PDF's as in Fig. 2. 

Figs. 4a and 4b refer to Au + Au at RHIC and to Pb + Pb at SPS respectively. We can clearly see that deviations 
from unity happen more strongly for RHIC collisions. 

The above situation has to be regarded as an extreme case, in which J/ip's are produced via parton fusion (with 
shadowing taken into account) in the initial state of the nucleus-nucleus collision and then just escape from the fireball, 
without further nuclear or comover suppression. 

We now consider the opposite extreme case, where all J/^'s are formed in the final state. Final state here means 
that the original J/^'s are totally destroyed by plasma screening and recriated at the final stage of the fireball life 
(at T — T c ~ 170 MeV) by coalescence and then just "emerge ready" during the fase transition, obeying the thermal 
distribution (42). This scenario is supported by the calculations performed in [3,17,42]. Assuming this production 
mechanism the ratio defined above is modified to: 



QCOM an Thermal I CCOM ' 



where Spgj^ = 0.25. When discussing thermal production a crucial aspect is how to normalize the distribution 
(42). The most straightforward procedure would be to normalize it in such a way that the total number of c and c 
quarks in the final state would match the number of charm quarks pairs produced initially and computed with the 
help of perturbative QCD. However, as we can see from Table I and also from the recent and comprehensive analysis 
performed in Ref. [18], already at the p + p — > cc + X level there is some uncertainty in the cross section coming 
from choices of the quark masses, renormalization scales and parton density parametrizations. A further source of 
uncertainty arises when we go from pp to AA collisions and introduce shadowing effects. Moreover, in a very recent 
calculation of thermal charm production [42], already taking into account the first PHENIX data [28], it was shown 
that, to explain the overall magnitude of the data, one needs to increase the NLO QCD cc yield by a factor 2.8. In 
view of this lack of precise knowledge, in our study we will normalize the integral of Eq. (42), Njj^ in two different 
ways: 

/) N W = 3 Nqq D ; II) N%* = 0.5 N$ D (49) 

where Nq^ d is given by Eq. (15), which, in turn, will be evaluated with two different nuclear parton distributions 

In Fig. 5 we plot the ratio Rm- The two lower (upper) pannels refer to Pb + Pb (Au + Au) data measured at SPS 
(RHIC). In the two left (right) pannels we use the Eskola et al. [38] (Hirai et al. [40]) nuclear parton distributions. 
Inside each of the pannels, the two upper curves are obtained with I) N'^ — 3 and the two lower curves with 

II) — 0.5 Nqq D . Finally, solid and dashed lines represent the maximal and the minimal contribution from the 

QGP component Eq. (39). Although f s = 10~ 2 it is easy to see that choosing a smaller value for the screening factor 
would move the solid lines to the dashed ones, since they contain already a very small number of "in-plasma" born 
J/V's. 

In Fig. 5 we see that there is already a rise of Rm at decreasing values of xp, but this rise has the shape of a 
plateau (as anticipated in the introduction), starting from xp ~ 0.01. If the screening is less effective, f s = 10~ 2 , then 
we find a deviation from the plateau at RHIC energies, as in Figs. 5a and 5b, but do not observe any visible effect 
at SPS energies, as in Figs. 5c and 5d. From these last figures we conclude that in-plasma production can hardly 
be detected at SPS since the xp distribution is always dominated by the thermal (plateau-like) contribution. On the 
other hand, it should be remarked that if a plateau is found experimentally, this is interesting in itself, being a strong 
evidence in favor of the statistical hadronization model, which implies plasma formation. 
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Now we come back to the prediction made in [4], namely that the initially produced J/i/>'s after the interaction 
either with hadronic comovers or with the plasma will disappear around xf — giving place to a dip in this region. 
We can split Eq. (47) into two pieces, the first one containing initial production and the second containing the plasma 
contribution. In order to include the suppression resulting from the final state interactions predicted in [4], we replace 
the first piece by a parametrization of Rab(%f, b — 0), Eq. (3) of [4], keeping our plasma contribution the same and 
construct a new version of (47): 

gCOM aN QGP 

r iv (*f) - r ab ( XF ) + „ F r jz^ (50) 

T A A(b = 0) dxp 

This quantity, plotted in Fig. 6 is a more realistic version of Rn (plotted in Fig. 4) . As in previous figures solid and 
dashed lines correspond to the nuclear parton densities of Refs. [38] and [40] respectively. Figs. 6a and 6b correspond 
to RHIC and SPS collisions respectively. Following the same convention employed in Fig. 4, in each pannel the 
two upper curves represent the largest QGP contribution, whereas the two lower curves represent the smallest QGP 
contribution. Comparing Figs. 4 and 6 we can see that the inclusion of FSI changes the results, decreasing the 
strength of the low xf enhancement. Nevertheless, the effect of in-plasma production remains quite visible. Notice 
that the lowest line crosses the origin of the horizontal axis (xj? ~ 0) at 0.16. This is still far above the point predicted 
in [4], of 0.055. 

We discuss now two more conservative scenarios, i.e., with less QGP production. First, we consider the lower part of 
Fig. 6a, amplifying the scale and allowing for smaller values of the screening factor f s (f s = 10~ 3 ). This amplification 
is shown in Fig. 7a). Inspite of the uncertainties in the calculations, our results in this figure, being still larger than 
the one found in [4] (the dash-dotted line), suggest that an enhancement in Riv in the low xf region might be seen 
if QGP would be formed. Finally, we compute again the QGP production rate, using the same inputs except for the 
initial temperature of the plasma, which we take to be To = 300 MeV. At these smaller temperatures the screening 
factor must be larger and we choose it to be f s = 0.1. This is shown in Fig. 7b) with the same conventions used in 
the previous figures. As it can be seen, our lowest line is still above the result of [4]. 



VII. CONCLUSIONS 



In this paper we have considered J/ip production in nucleus-nucleus collisions by three different mechanisms: direct 
(primordial parton fusion) ii) thermal (statistical coalescence at hadronization) and iii) QGP (in-plasma parton fusion). 
As already discussed in [3] i) and ii) are just two extreme cases. In realistic simulations, for a fixed number of cc pairs 
and of J/ip's, these last are gradually destroyed, giving origin to the "regenerated" component ii). Component iii) is 
reconsidered here after being forgotten for long time. Our main point was that, even being small, contribution iii) is 
very strongly peaked around xp — and can thus become visible if enough plasma is formed. 

With this in mind we defined the ratio (46) and computed it assuming mechanism i) plus a small component of iii) 
(ratio Rn) and assuming mechanism ii) plus a small component of iii) (ratio Rm). In the first case it is completely 
flat without plasma contribution. A QGP contribution creates a big bump at xf < 0.01. In the second case, without 
plasma contribution, we observe a step structure, with a plateau at xp < 0.01. Switching on the plasma contribution 
creates a steeply falling curve. We believe that any of these features can be measured and will be very interesting for 
charm physics at RHIC. 

In [4] it was suggested that the xf distribution to be measued at RHIC is flat and has a dip in the region xf — 0. 
Our calculations indicate that the inclusion of the "in-plasma born" J/ip's will fill this dip. 

Of course, there are several improvements that one could do in the in-plasma production, as, for example, the use 
of a more sophisticated hydrodynamical expansion. Moreover, final state interactions of this component should be 
considered. Work along these lines is already in progress. Our findings in this work encourage us to pursue this 
program. 
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QCD 


m c = 1.2 


GeV 


m c = 1.5 


GeV 




Eskola et al. [38] 


Hirai et al. [40] 


Eskola et al. [38] 


Hirai et al. [40] 


Nqcd (without shadowing) 


6.52 


4.65 


2.65 


2.07 


(with shadowing) 


6.32 


4.36 


2.61 


1.97 


Nqq D (without shadowing) 


0.096 


0.063 


0.019 


0.013 


(with shadowing) 


U.Uoo 


0.058 


0.018 


0.012 


TABLE I. Number of J ftp's 


and cc pairs produced in Au 


— Au collisions at RHIC from QCD calculations for different values 


of couplings and charm quark 


masses. 








DPP 


Tfl c 


1.2 GeV 


m c = 


1.5 GeV 












Levai el al. [6] 




- (3.7) 




- (1-1) 


n qgp Rafelski el al. [32] 






- (15) 




Miiller el al. [7] 






17 - 




This work 


120 (60) 


39(19.5) 


22 (11) 


7.6(3.8) 


N QGP This work 


2.70 


0.84 


0.49 


0.16 



TABLE II. Number of "in-plasma" produced cc pairs and J ftp's in RHIC collisions. Numbers inside parenthesis are obtained 
with k — 1 and the others with k = 2. The initial temperature is 550 MeV in all cases. 
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a) b) 




FIG. 1. J/ip momentum distributions: direct QCD (13) (dash-dotted line), in-plasma (39) (solid line) and thermal production 
(42) (dashed line). All curves have the same normalization (see text), a) RHIC collisions; b) SPS collisions. 
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FIG. 2. Ratio Ri(xf): effect of shadowing in QCD direct production. The numerators were obtained with (13) and (12) 
and the nuclear parton distribution functions were taken from [38] (solid lines) and from [40] (dashed lines). 
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FIG. 5. Ratio Riii(xf) calculated with (48); labels I and II represent the normalizations defined in (49); solid (dashed) 
lines are obtained with the largest (smallest) in-plasma contribution; a) and b) are for RHIC collisions and c) and d) for SPS 
collisions; in the two left (right) pannels we use the Eskola et al. [38] (Hirai et al. [40]) nuclear parton distributions. 
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FIG. 7. Ratio Riv(xf) computed with (50) for RHIC collisions; solid lines: parton densities of Ref. [38]; dashed lines: parton 
densities of Ref. [40] ; dash-dotted line represents the (parameterized) result for central collisions obtained by Hiifner et al. [4] . 
7a) To = 550 MeV; 7b) T = 300 MeV. 
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